TCGA BRCA


TCGA-BRCA data preparation

Download TCGA-BRAC data

library(TCGAbiolinks)

GDCquery link

# query <- GDCquery(project = "TCGA-BRCA",
#                   data.category = "Transcriptome Profiling",
#                   data.type = "Gene Expression Quantification")
# 
query <- GDCquery(project = "TCGA-BRCA",
                  data.category = "Transcriptome Profiling",
                  data.type = "Gene Expression Quantification", 
                  sample.type = c("Primary Tumor","Solid Tissue Normal"))
df = query$results[[1]]
df$sample_type %>% table() %>% data.frame()
# Download and prepare data
GDCdownload(query)
data <- GDCprepare(query)
library(DESeq2)
# 가능한 모든 assays 이름을 확인
assayNames(data)
# unstranded: Strand 비특이적 발현 데이터
# stranded_first: Strand 특이적 (첫 번째 strand) 발현 데이터
# stranded_second: Strand 특이적 (두 번째 strand) 발현 데이터
# tpm_unstrand: TPM (Transcripts Per Million) 방식으로 정규화된 strand 비특이적 발현 데이터
# fpkm_unstrand: FPKM (Fragments Per Kilobase of transcript per Million mapped reads) 방식으로 정규화된 strand 비특이적 발현 데이터
# fpkm_uq_unstrand: FPKM-UQ (Upper Quartile normalized FPKM) 방식으로 정규화된 strand 비특이적 발현 데이터

Select the raw reads

# 'unstranded' 데이터를 사용하여 발현 수치 추출
counts <- assay(data, "unstranded")
# ENSG to Symbols
library(org.Hs.eg.db)

# 유전자 ID에서 버전 제거
ensembl_ids <- gsub("\\..*", "", rownames(counts))

# 중복된 ID 처리
unique_ids <- unique(ensembl_ids)

# 중복 제거된 ID에 대해 유전자 심볼 매핑
gene_symbols <- mapIds(org.Hs.eg.db,
                       keys = unique_ids,
                       column = "SYMBOL",
                       keytype = "ENSEMBL",
                       multiVals = "first")

# 원래 데이터에 매핑된 유전자 심볼 할당
# 여기서 중복 제거된 목록을 이용해 각각의 원본 ID에 매핑된 심볼 할당
symbol_names <- gene_symbols[ensembl_ids]
rownames(counts) <- symbol_names
# Metadata 
filtered_colData <- colData(data)[colnames(counts),]
TCGA_BRCA_countsMeta = list(counts = counts,
                            meta = filtered_colData)
TCGA_BRCA_countsMeta %>% saveRDS(paste0(dir,"TCGA_BRCA_countsMeta.rds"))

TCGA BRCA 프로젝트의 메타데이터

제공된 열 이름 목록은 TCGA BRCA 프로젝트의 메타데이터에서 사용된 모든 열을 나타냅니다. 각 열의 용도는 다양하며, 환자와 샘플에 관련된 상세한 정보를 포함합니다. 이 정보는 연구자가 종양의 특성, 환자의 임상 경로, 생물학적 데이터와의 관계를 분석하는 데 사용됩니다. 아래는 각 열 이름의 의미와 사용 용도를 설명합니다:

  1. barcode: TCGA 프로젝트 내 각 샘플에 대한 고유 식별 코드입니다.
  2. patient: 각 환자에게 할당된 고유 식별자입니다.
  3. sample: 샘플의 고유 식별자입니다.
  4. shortLetterCode: 샘플의 간단한 코드를 나타냅니다.
  5. definition: 샘플 유형의 정의를 설명합니다.
  6. sample_submitter_id: 제출된 샘플의 ID입니다.
  7. sample_type_id: 샘플 유형의 ID입니다.
  8. tumor_descriptor: 종양의 설명을 제공합니다.
  9. sample_id: 샘플 ID입니다.
  10. sample_type: 샘플의 유형을 나타냅니다 (예: Primary Tumor).
  11. composition: 샘플 구성을 설명합니다.
  12. days_to_collection: 샘플 수집까지 걸린 일수입니다.
  13. state: 샘플의 상태를 나타냅니다.
  14. initial_weight: 샘플의 초기 무게입니다.
  15. preservation_method: 샘플 보존 방법입니다.
  16. pathology_report_uuid: 병리 보고서의 고유 식별자입니다.
  17. submitter_id: 제출자 ID입니다.
  18. oct_embedded: OCT(광학 단층 촬영)에 포함되었는지 여부입니다.
  19. specimen_type: 표본 유형입니다.
  20. is_ffpe: FFPE(형식고정 파라핀 포매) 처리 여부입니다.
  21. tissue_type: 조직 유형입니다.
  22. synchronous_malignancy: 동시 악성 종양의 유무입니다.
  23. ajcc_pathologic_stage: AJCC 병리학적 단계입니다.
  24. days_to_diagnosis: 진단까지 걸린 일수입니다.
  25. treatments: 처리 내용입니다.
  26. last_known_disease_status: 마지막으로 알려진 질병 상태입니다.
  27. tissue_or_organ_of_origin: 기원 조직 또는 기관입니다.
  28. days_to_last_follow_up: 마지막 추적 관찰까지 걸린 일수입니다.
  29. age_at_diagnosis: 진단 당시 나이입니다.
  30. primary_diagnosis: 주 진단입니다.
  31. prior_malignancy: 이전 악성 종양의 유무입니다.
  32. year_of_diagnosis: 진단 연도입니다.
  33. prior_treatment: 이전 치료 내용입니다.
  34. ajcc_staging_system_edition: AJCC 병기 시스템 에디션입니다.
  35. ajcc_pathologic_t: AJCC 병리학적 T 등급입니다.
  36. morphology: 형태학입니다.
  37. ajcc_pathologic_n: AJCC 병리학적 N 등급입니다.
  38. ajcc_pathologic_m: AJCC 병리학적 M 등급입니다.
  39. classification_of_tumor: 종양 분류입니다.
  40. diagnosis_id: 진단 ID입니다.
  41. icd_10_code: 국제 질병 분류 코드(ICD-10)입니다.
  42. site_of_resection_or_biopsy: 절제 또는 생검 부위입니다.
  43. tumor_grade: 종양 등급입니다.
  44. progression_or_recurrence: 진행 또는 재발 여부입니다.
  45. alcohol_history: 음주력입니다.
  46. exposure_id: 노출 ID입니다.
  47. race: 인종입니다.
  48. gender: 성별입니다.
  49. ethnicity: 민족입니다.
  50. vital_status: 생존 상태입니다.
  51. age_at_index: 인덱스 시 나이입니다.
  52. days_to_birth: 출생까지 걸린 일수입니다.
  53. year_of_birth: 출생 연도입니다.
  54. demographic_id: 인구통계학적 ID입니다.
  55. days_to_death: 사망까지 걸린 일수입니다.
  56. year_of_death: 사망 연도입니다.
  57. bcr_patient_barcode: BCR 환자 바코드입니다.
  58. primary_site: 주요 발생 부위입니다.
  59. project_id: 프로젝트 ID입니다.
  60. disease_type: 질병 유형입니다.
  61. name: 이름입니다.
  62. releasable: 공개 가능 여부입니다.
  63. released: 공개 여부입니다.
  64. days_to_sample_procurement: 샘플 조달까지 걸린 일수입니다.
  65. paper_patient: 논문에 사용된 환자 데이터입니다.
  66. paper_Tumor.Type: 논문에 기술된 종양 유형입니다.
  67. paper_Included_in_previous_marker_papers: 이전 마커 논문에 포함되었는지 여부입니다.
  68. paper_vital_status: 논문에 기술된 생존 상태입니다.
  69. paper_days_to_birth: 논문에 기술된 출생까지의 일수입니다.
  70. paper_days_to_death: 논문에 기술된 사망까지의 일수입니다.
  71. paper_days_to_last_followup: 논문에 기술된 마지막 추적 관찰까지의 일수입니다.
  72. paper_age_at_initial_pathologic_diagnosis: 논문에 기술된 최초 병리학적 진단 시 나이입니다.
  73. paper_pathologic_stage: 논문에 기술된 병리학적 단계입니다.
  74. paper_Tumor_Grade: 논문에 기술된 종양 등급입니다.
  75. paper_BRCA_Pathology: 논문에 기술된 BRCA 병리학입니다.
  76. paper_BRCA_Subtype_PAM50: 논문에 기술된 BRCA 하위 유형 PAM50입니다.
  77. paper_MSI_status: 논문에 기술된 MSI 상태입니다.
  78. paper_HPV_Status: 논문에 기술된 인유두종바이러스(HPV) 상태입니다.
  79. paper_tobacco_smoking_history: 논문에 기술된 흡연력입니다.
  80. paper_CNV Clusters: 논문에 기술된 CNV 클러스터입니다.
  81. paper_Mutation Clusters: 논문에 기술된 돌연변이 클러스터입니다.
  82. paper_DNA.Methylation Clusters: 논문에 기술된 DNA 메틸화 클러스터입니다.
  83. paper_mRNA Clusters: 논문에 기술된 mRNA 클러스터입니다.
  84. paper_miRNA Clusters: 논문에 기술된 miRNA 클러스터입니다.
  85. paper_lncRNA Clusters: 논문에 기술된 lncRNA 클러스터입니다.
  86. paper_Protein Clusters: 논문에 기술된 단백질 클러스터입니다.
  87. paper_PARADIGM Clusters: 논문에 기술된 PARADIGM 클러스터입니다.
  88. paper_Pan-Gyn Clusters: 논문에 기술된 Pan-Gyn 클러스터입니다.

TCGA-BRCA Analysis

TCGA_BRCA_countsMeta = readRDS(paste0(dir,"TCGA_BRCA_countsMeta.rds")) 

Create DESeq object

count.raw = TCGA_BRCA_countsMeta$counts

# Remove the duplicated gene names 
rs = rownames(count.raw)[!is.na(rownames(count.raw))]
count.mtx = count.raw[rs,]

# Remove genes with no expression across samples 
count.mtx = count.mtx[rowSums(count.mtx) !=0,]
meta = TCGA_BRCA_countsMeta$meta
meta = meta[colnames(count.mtx),]
se <- SummarizedExperiment(as.matrix(count.mtx), 
                           colData=meta)
dds <- DESeqDataSet(se, ~ 1)
vsd <- vst(dds, blind=FALSE)

pcaData <- DESeq2::plotPCA(vsd, intgroup = "sample_type", returnData = TRUE)
pcaData$group = dds$group
PCA_var=attr(pcaData, "percentVar")

PCA plot

ggplot(pcaData, aes(x = PC1, y = PC2, fill = sample_type)) +
  geom_point(size = 2, alpha = 0.8, shape = 21, color = "black", stroke = 0.2)  +
  # ggrepel::geom_text_repel(aes(label=name), 
  #                          color="grey6", size=3, hjust= -0.3, vjust=-0.3) +
  labs(x = paste("PC1: ", round(100 * PCA_var[1]), "% variance"),
       y = paste("PC2: ", round(100 * PCA_var[2]), "% variance")) +
  theme_bw() +
  theme(legend.title = element_blank()) +
  ggtitle("PCA") +
  labs(caption = " ")

Differentially Expressed Genes (DEG) analysis

DEG strategy

Differentially Expressed Genes (DEG) between Tumor tissue and Normal tissue

Tumor tissue/Normal tissue

n = 1224
Tumor n = 1111
Normal n = 113




Create DEG object

# Generate info table
info <- data.frame(matrix(nrow = ncol(count.mtx), ncol = 2))
colnames(info) <- c('sample', 'cond')
info$sample <- colnames(count.mtx)
info$cond <- dds$sample_type
info$cond <- factor(info$cond, 
                    levels = c("Solid Tissue Normal","Primary Tumor")) # CTL going first
# levels(info$cond)

# DESeq
dds <- DESeqDataSetFromMatrix(count.mtx, info, ~ cond)
dds <- DESeq(dds) 
# dds %>% saveRDS(paste0(dir,"TCGA_BRCA_countsMeta.dds.rds"))
res <- results(dds)
res <- data.frame(res)
dds = readRDS(paste0(dir,"TCGA_BRCA_countsMeta.dds.rds"))
res <- results(dds)
res <- data.frame(res)
# Add DEG information 
fc = 1.5
pval = 0.05

# res = res %>% mutate(DE=ifelse(log2FoldChange >= log2(fc) & pvalue < pval, 'UP',
#                                ifelse(log2FoldChange <= -log2(fc) & pvalue < pval, 'DN','no_sig')))

res = res %>% mutate(DE=ifelse(log2FoldChange >= log2(fc) & padj < pval, 'UP',
                               ifelse(log2FoldChange <= -log2(fc) & padj < pval, 'DN','no_sig')))
res = na.omit(res)

DEG table

res %>% DT::datatable()

Volcanoplot

res$DE = factor(res$DE, levels = c('UP','DN','no_sig'))
res %>% 
  ggplot(aes(log2FoldChange, -log10(padj), color=DE)) + 
  geom_point(size=1, alpha=0.5) + 
  scale_color_manual(values = c("red3","royalblue3","grey"), guide = FALSE) +
  theme_classic() +
  geom_vline(xintercept = c(-log2(fc),log2(fc)), color='grey') +
  geom_hline(yintercept = -log10(0.05),color='grey') +
  guides(colour = guide_legend(override.aes = list(size=5))) +
  ggtitle(paste0(levels(dds$cond)[2], " / ", levels(dds$cond)[1] )) +
  ggeasy::easy_center_title() ## to center title

Volcanoplot with Number of DEGs

t= paste0(levels(dds$cond)[2], " / ", levels(dds$cond)[1] )
up = nrow(res[res$DE == "UP", ])
dn = nrow(res[res$DE == "DN", ])
res %>% ggplot(aes(log2FoldChange, -log10(padj), color=DE)) + 
  geom_point(size=0.5, shape=19, alpha=0.7) +
  geom_vline(xintercept = c(-log2(fc), log2(fc)), size=0.1, color="grey") +
  geom_hline(yintercept = -log10(0.05), size=0.1, color="grey") +
  scale_color_manual(values = c("red3","royalblue3","grey"), guide = FALSE) +
  theme_bw() +
  annotate("text", x = Inf, y = Inf, label = paste0("UP: ", up), 
           hjust = 1.1, vjust = 2, size = 5, color = "red") +
  annotate("text", x = -Inf, y = Inf, label = paste0("DN: ", dn), 
           hjust = -0.1, vjust = 2, size = 5, color = "royalblue") +
  theme_bw() + ggtitle(t)

Number of DEGs

res %>% filter(DE != "no_sig") %>% 
  ggplot(aes(DE, fill=DE)) + geom_bar(color="black", size=0.2) +
  geom_text(stat = 'count', aes(label = ..count..), vjust = -0.1, size= 4, color=c("salmon","royalblue")) +
  scale_fill_manual(values = c("salmon", "royalblue"), guide=F) +
  theme_bw()

GSEA function

library(clusterProfiler)
hallmark <- msigdbr::msigdbr(species = "Homo sapiens", category = "H") %>% 
  dplyr::select(gs_name, gene_symbol)

perform_GSEA <- function(res, ref, pvalueCutoff = 1) {
  ranking <- function(res) {
    df <- res$log2FoldChange
    names(df) <- rownames(res)
    df <- sort(df, decreasing = TRUE)
    return(df)
  }
  
  ranked.res <- ranking(res)
  
  x <- clusterProfiler::GSEA(geneList = ranked.res,
                             TERM2GENE = ref,
                             pvalueCutoff = pvalueCutoff,
                             pAdjustMethod = "BH",
                             verbose = TRUE,
                             seed = TRUE)
  
  result <- x@result %>% arrange(desc(NES))
  result <- result[, c('NES', 'pvalue', 'p.adjust', 'core_enrichment', 'ID')]
  return(result)
}

# Application 
gsea.res = perform_GSEA(res = res, ref = hallmark) 
filtered_gsea = gsea.res %>% mutate(sig= ifelse(pvalue <= 0.05,"p value <= 0.05", "p value > 0.05"))

# Modified GSEA NES plot 
gsea_nes_plot =function(gsea.res, title, fontsize.x = 5, fontsize.y = 6){
  gsea.res %>% ggplot(aes(reorder(ID, NES), NES)) +
    geom_col(aes(fill=sig), color="grey1", size=0.2) +
    coord_flip() +
    labs(x="Pathway", y="Normalized Enrichment Score",
         title= "GSEA") + 
    theme_classic() +
    # scale_fill_gradient(low = '#FF0000', high = '#E5E7E9') +
    scale_fill_manual(values = c("#FF0000","grey88")) +
    theme(axis.text.x= element_text(size=fontsize.x, face = 'bold'),
          axis.text.y= element_text(size=fontsize.y, face = 'bold'), 
          axis.title =element_text(size=10)) +ggtitle(title)
}

GSEA NES plot

t= paste0("DEGs from ",levels(dds$cond)[2], " / ", levels(dds$cond)[1] )
gsea_nes_plot(filtered_gsea, title = t, fontsize.x = 8, fontsize.y = 8)

GSEA table

gsea.res %>% DT::datatable()